The present invention relates to a method of simulating a semiconductor device.
The Semiconductor device simulation techniques have been known and disclosed for examples in "Process Device Simulation Technique" Industrial Book, pp. 99-104, 1990, R. Thoma et al. "Hydrodynamic Equations For Semiconductors with Non-parabolic Band Structure", IEEE Transactions on Electron devices, vol. 38, No. 6, pp. 1343-1353, June 1991, and in "A Study on Relaxation Time Application for Semiconductor Device Analysis" SSD 84-68, Vol. 84, No. 180, pp. 31-36, 1984, as well as in K. W. Chai et al. "Hydrodynamic Simulation Of Electron Heating In Conventional And Lightly-Doped-Drain MOSFETs With Application To Substrate Current Calculation" International Journal of Numerical Modelling: Electronic Networks, Devices and Fields, vol. 5, pp. 53-66, 1992 and further in R. K. Cook, "Numerical Simulation Of Hot-Carrier Transport In Silicon Bipolar Transistor", IEEE Transaction On Electron Devices, Vol. 30, No. 9, pp. 1103-1109, 1983.
An outline of the general device simulation will be described as follows. In numerical analysis to semiconductor devices, a drift-diffusion model and an energy transport model may be used. In the drift-diffusion model, carriers, for example, electrons and holes are approximated as fluids. The energy transport model uses approximation in higher order. The device simulation is made in the drift-diffusion model in a steady state, wherein the following equation of conservation of charges, electron current continuous equation, and hole current continuous equation are set. EQU div D=.rho. (equation of conservation of charges) (1) EQU D=.epsilon. E (2) EQU E=-grad.upsilon. (3) EQU .rho.=q(p-n+N.sub.D -N.sub.A) (4)
D:electric displacement PA0 .rho.: density of charge PA0 E: electric field PA0 .epsilon.: dielectric constant PA0 q: magnitude of electronic charge PA0 p: hole concentration PA0 n: electron concentration PA0 N.sub.D : donor concentration PA0 N.sub.A : acceptor concentration EQU div Jn=q(R-G) (electron current continuous equation) (5) EQU div Jp=-q(R-G) (hole current continuous equation) (6) PA0 Jn: electron current PA0 Jp: hole current PA0 R: carrier recombination term PA0 G: carrier generation term EQU Jn=q n .mu.n E+q Dn grad (n) (7) EQU Jp=q p .mu.p E-q Dp grad (p) (8) PA0 .mu.n: electron mobility PA0 .mu.p: hole mobility PA0 Dn: electron diffusion coefficient PA0 Dp: hole diffusion coefficient EQU Dn=.mu.n {(k.sub.B T)/q} (9) EQU Dp=.mu.p {(k.sub.B T)/q} (10) PA0 k.sub.B : Boltzmann's constant PA0 T: lattice temperature PA0 T*.sub.n : electron temperature PA0 T*.sub.p : hole temperature PA0 .tau..sub.in : electron momentum relaxation time PA0 .tau..sub.ip : hole momentum relaxation time EQU .tau.*.sub.in =(1/3)m*.sub.n (M.sub.n.sup.-1) .tau..sub.in (19) EQU .tau.*.sub.ip =(1/3)m*.sub.p (M.sub.p .sup.-1) .tau..sub.ip (20) PA0 m*.sub.n : electron effective mass PA0 m*.sub.p : hole effective mass PA0 M.sub.n.sup.-1 : electron inversion effective mass tensor PA0 M.sub.p.sup.-1 : hole inversion effective mass tensor PA0 (): average operation in k space EQU div Sn=-Jn grad.upsilon.-(3/2)k.sub.B n {(T*.sub.n -T.sub.neq)/.tau.*wn}(21) PA0 (electron energy conservation) EQU div Sp=Jp grad.upsilon.-(3/2)k.sub.B p {(T*.sub.p -T.sub.pcq)/.tau.*wp}(22) PA0 (hole energy conservation) PA0 Sn: electron energy flow density PA0 Sp: hole energy flow density PA0 T.sub.neq : electron equilibrium temperature PA0 T.sub.peq : hole equilibrium temperature EQU .tau.*.sub.wn =(3/2)k.sub.B (T*.sub.n -T.sub.eq){T.sub.wn ( (.epsilon..sub.n) (.epsilon..sub.neq) )} (23) EQU .tau.*.sub.wp =(3/2)k.sub.B (T*.sub.p -T.sub.eq){T.sub.wp ( (.epsilon..sub.p) (.epsilon..sub.peq) )} (24) PA0 (.epsilon..sub.n): mean electron energy PA0 (.epsilon..sub.p): mean hole energy PA0 (.epsilon..sub.neq): electron equilibrium energy PA0 (.epsilon..sub.peq): hole equilibrium energy PA0 .tau.*.sub.wn : electron energy relaxation time PA0 .tau.*.sub.wp : hole energy relaxation time EQU Sn=-5/2(k.sub.B T*.sub.n /q) (.tau.*.sub.sn /.tau.*.sub.in) {Jn+(q/m*.sub.n) .tau..sub.in n grad (k.sub.B T*.sub.n)} (25) EQU Sp=-5/2(k.sub.B T*.sub.p /q) (.tau.*.sub.sp /.tau.*.sub.ip) {Jp+(q/m*.sub.p) .tau..sub.ip p grad (k.sub.B T*.sub.p)} (26) EQU .tau.*.sub.sn =(1/3(M.sub.n.sup.-1 .epsilon..sub.n +v.sub.n v.sub.n) )/(5/6(v.sub.n.sup.2) ) .tau..sub.sn (27) EQU .tau.*.sub.sp =(1/3(M.sub.p.sup.-1 .epsilon..sub.p +v.sub.p v.sub.p) )/(5/6(v.sub.p.sup.2) ) .tau..sub.sp (28) PA0 .tau..sub.sn : relaxation time for electron energy flow density Sn PA0 .tau..sub.sp : relaxation time for hole energy flow density Sp PA0 Vn: electron velocity PA0 Vp : hole velocity
In the above equations, variables to be solved are the potential .upsilon., the electron concentration n, and hole concentration p.
For the energy transport model in the steady state, the following equation is set wherein the equation of the drift-diffusion model is added with the carriers energy conservation equations. EQU div D=.rho. (equation of conservation of charges) (11) EQU D=.epsilon. E (12) EQU E=-grad.upsilon. (13) EQU .rho.=q(p-n+N.sub.D -N.sub.A) (14) EQU div Jn=q(R-G) (electron current continuous equation) (15) EQU div Jp=-q(R-G) (hole current continuous equation) (16) EQU Jn=q n.mu.n E+.mu.n (.tau..sub.in /.tau.*.sub.in) grad (n k.sub.B T*.sub.n) (17) EQU Jp=q p.mu.p E-.mu.p (.tau..sub.ip /.tau.*.sub.ip) grad (p k.sub.B T*.sub.p) (18)
In the above equations in energy transport equation model, variables to be solved are potential .upsilon., electron concentration n, hole concentration p, electron temperature T*.sub.n, and hole temperature T*.sub.p. Whereas the mark * has been put to distinguish the above temperatures from the temperatures as defined in the thermodynamics, such mark will be omitted hereinafter. EQU Tn=2/(3 k.sub.B) (.epsilon..sub.n) (29) EQU Tp=2/(3 k.sub.B) (.epsilon..sub.p) (30)
Biases are sequentially renewed by setting designated plural applied devices as boundary conditions in order to calculate the five equations of equation of conservation of charges, electron current continuous equation, hole current continuous equation, electron energy conservation equation and hole energy conservation equation.
Since those equations are non-linear equations, the solutions are found by Newton method as follows.
If the following equation (31) has been given for variable "x" and addition of a variable .delta.x0 to x0 provides the solution, then the following equation (32) is given. EQU F(x)=0 (31) EQU F(x0+.delta.x0)=0 (32)
The differential coefficient of F(x) is defined to be F'(x0) and Taylor expansion is made in the first order to F(x0+.delta.x0) for .delta.x0 to give the following equations. EQU F(x0+.delta.x0)=F(x)+F'(x0).delta.x0=0 (33) EQU .delta.x0=-F(x0)/F'(x0) (34)
The following equation (35) is set for similar calculation for x1. EQU X1=X0+.delta.x0 (35)
The above calculation will be repeated so that if .delta.xi is smaller than convergence condition .epsilon., this means convergence. xi is considered to be the solution of the equation (31).
FIG. 1 is a flow chart illustrative of the above calculations. FIG. 2, is a graph illustrative of a relationship of F(x) and x. x approaches from the initial value x0 to the solution xi. Tangential lines with the function curve y=F(x) are set to find intersection points of the tangential lines and the x-axis. The intersection point is set to be the next value of x. If the initial value x0 is close to the solution, then the number of the required interating calculations is small and it takes a short time to obtain the solution. This means that the initial value x0 closer to the solution is better initial value.
Whereas in the above descriptions to the Newton method, the single variable is dealt with, in the device simulation, meshes are formed over entire analysis regions so that equations are set for variables on the meshes. FIG. 3 is a view illustrative of a mesh for analysis. Each of the five parameters, for example, potential, electron concentration, hole concentration, electron temperature and hole temperature represents variables, the number of which corresponds to the number of the points N of the mesh. This means that 5N simultaneous equations are solved. The above charge conservation equation, electron current continuous equation, hole current continuous equation, electron energy conservation equation and hole energy conservation equation are expressed as follows. EQU F.upsilon.(.upsilon., n, p, Tn, Tp)=0 (charge conservation equation) (36) EQU Fn(.upsilon., n, p, Tn, Tp)=0 (electron current continuous equation) (37) EQU Fp(.upsilon., p, Tn, Tp)=0 (hole current continuous equation) (38) EQU FTn(.upsilon., n, p, Tn, Tp)=0 (electron energy conservation equation) (39) EQU FTp(.upsilon., n, p, Tn, Tp)=0 (hole energy conservation equation) (40)
where .upsilon. is the potential, n is the electron concentration, p is the hole concentration, Tn is the electron temperature, and Tp is the hole temperature, and each represents N variables.
In the above cause, there may be used a coupled method and a de-coupled method. In accordance with the coupled method, the above five simultaneous equations are concurrently solved for all of the variables. By contrast, in accordance with the de-coupled method, the above five equations are separately solved for individual variables.
FIG. 4 is a flow chart illustrative of the procedure of the coupled method. FIG. 5 is a flow chart illustrative of the procedure of the decoupled method. F.upsilon.n' represents the following partial differential. EQU F.upsilon.n'=(.differential.F.upsilon./.differential.n) (41)
In accordance with the coupled method, the above five simultaneous equations are concurrently solved for all of the variables.
On the other hands, in accordance with the de-coupled method, all variables except for the currently attracted variable are fixed to solve the equations. If, for example, the electron energy conservation equation is solved by setting electron temperature as variable but fixing the potential, the electron concentration, hole concentration and hole temperature.
In each iteration, the matrix is calculated. In the coupled method, for the single iteration, a single matrix of N.times.N is solved. By contrast, in the de-coupled method, five matrixes of N.times.N are solved.
The coupled method is capable of finding the solution through a small number of iterations. The result of the calculation by the coupled method is largely dependent upon the initial value. Notwithstanding, if the initial value is not proper, then it might be difficult to obtain the convergence.
The de-coupled method is not so dependent upon the initial value, but requires a large number of the iterations of calculations. It takes a shorter time to calculate in a single iterations by the de-coupled method as compared to when the calculation is made by the coupled method. The decoupled method requires a large number of the iterations of the calculations as compared to when the calculation is made by the coupled method.
In the majority case, the required calculation time for obtaining the solution is shorter when the calculation is made by the coupled method as compared to when the calculation is made by the de-coupled method. If the better initial value can be given, it can take a shorter time to analyze the semiconductor device through the coupled method. This means that it is important to set the better initial value for calculation through the coupled method.
A control volume method may be used for discretization from the equations into the analysis mesh. FIG. 6 is a diagram illustrative of the polygon which comprises a plurality of triangle meshes which are represented by real lines. Bisectors of individual mesh edges of the triangle meshes form another polygon which represents the control volume. Vertexes of the control volume polygon correspond to circumcenters of the triangle meshes, namely centers of the circumscribing circles.
In the control volume method, the flow of the physical quantities such as current can be expressed by the product of the density of the physical quantity flow over the mesh edge such as current density and the length of the side of the control volume polygon. In the two dimensional case, in place of the length of the side of the control volume polygon, cross section is used.
The following description will focus on the prior art for estimation of the initial value of the energy transport model analysis. For the energy transport model analysis, it is preferable to set initial values of an electron temperature and a hole temperature to be closer to the solutions.
In a first conventional estimation method, the thermal equilibrium temperatures are set as initial values of the carrier temperatures and may expressed by the following equations. FIG. 7 is a flow chart illustrative of the procedures of the first conventional estimation method. EQU Tnk=Tneq (k=1-N) (42) EQU Tpk=Tpeq (k=1-N) (43)
In the energy transport model analysis, as a bias is applied to electrode of the semiconductor device, the carrier temperatures are increased to become different from the thermal equilibrium temperatures, for which reason the initial value is likely to be set far from the solution under the bias applied condition, resulting in a deterioration of the convergence.
In a second conventional estimation method, the analysis is approximately made to obtain the initial value. The procedures will be described with reference to FIG. 8 which is a flow chart illustrative of the procedures of the second conventional estimation method. In a first step 1201, a bias condition is set. In a second step 1202, the drift-diffusion model is solved under the bias condition to obtain a potential, an electron concentration and a hole concentration. In a third step 1203, the obtained potential and electron and hole concentrations are used to solve the electron energy conservation equation and the hole energy conservation equation thereby to obtain an electron temperature and a hole temperature. In a fourth step 1204, the obtained electron temperature and hole temperature are defined to be carrier temperature initial values.
In the above conventional method, the energy conservation equation is solved one time, for which reason it is possible to obtain the initial value relatively close to the solution. Since, however, calculation is made by use of the matrix to solve the energy conservation equation, it takes a long time to make the required calculation.
In a third conventional estimation method, the solution having been obtained by analysis under the previous bias condition is set to be the initial value. FIG. 9 is a flow chart illustrative of the procedures of the first conventional estimation method. In a first step 1301, for initial bias conditions, thermal equilibrium temperatures are set to be initial values of the carrier temperatures to calculate the energy transport model as the transitional analysis to the applied bias. In a second step 1302, for the second and later bias conditions, the solutions having been obtained by the analysis to the previous bias conditions are set to be the carrier temperature initial values.
In the third conventional method, no calculation to the matrix is made for setting the initial value, for which reason it takes a short time to make the required calculation for the device analysis. If the applied biases have been prepared precisely, it is possible to set the desired initial value which is close to the solution. If, however, the variation in the applied bias is large, the initial value may be far from the solution, resulting in the deterioration of the convergence.
As described above, the first conventional method of estimation of the initial value for energy transport model analysis is engaged with the first problem with deterioration of the convergence under the highly applied bias conditions.
The second conventional method of estimation of the initial value for energy transport model analysis is engaged with the second problem that it takes a long time to make the required calculation for the required analysis.
The third conventional method of estimation of the initial value for energy transport model analysis is engaged with the third problem with deterioration of the convergence under the conditions of the large variation in applied bias.
In the above circumstances, it had been required to develop a novel method of simulating a semiconductor device free from the above problems.